# -*- coding: utf-8 -*-
"""
Created on Mon Jul  5 15:25:56 2021

@author: os4875st
"""
import numpy as np
import matplotlib.pyplot as plt
import cv2
from skimage import io, exposure
from skimage.util.dtype import img_as_float32
import fun_figs as ff
import tifffile
import time
import os, sys
import math
import copy

#Personal files
import file_processing as fp
import tiff_image_manipulation as ti
import functions_nd2_handling
import functions_pillar_free_wave_generation as ap

def align_two_stacks_manually(s1, s2, save_as_uint16, dx=0, dy=0):
    print('\tAligning stacks manually')
    s2_new = np.zeros_like(s2, dtype=s2.dtype)
    n, h, w = s2.shape
    #x-direction
    if dx > 0:
        x1 = 0+dx
        x2 = w
        d = x2-x1
        x1_org = 0
        x2_org = w-dx
    elif dx <= 0:
        x1 = 0
        x2 = w+dx
        d = x2-x1
        x1_org = w-d
        x2_org = w    

    #y-direction
    if dy > 0:
        y1 = 0+dy
        y2 = h
        d = y2-y1
        y1_org = 0
        y2_org = h-dy
    elif dy <= 0:
        y1 = 0
        y2 = h+dy
        d = y2-y1
        y1_org = h-d
        y2_org = h  

    s2_new[:,y1:y2,x1:x2] = s2[:,y1_org:y2_org,x1_org:x2_org]
    ff.show_img(s2_new[0], title_text = 'After alignment\n(Right image)', cmap='gray', sz=5, fsz = 15)
    return s2_new 

def align_two_stacks(s1, s2, save_as_uint16):
    """
    Translational alignment between two image stacks
    See https://learnopencv.com/image-alignment-ecc-in-opencv-c-python/ to learn more
    Parameters
    ----------
    s1 : (numpy array)
        Dimensions should be ordered ZYX
    s2 : (numpy array)
        Dimensions should be ordered ZYX

    Returns
    -------
    s2_corrected : (numpy array, uint8)
        s2 aligned to s1

    Updates: try uint8 instead of float32?
            Reduce termination_eps if the code is slow, see https://stackoverflow.com/questions/45997891/cv2-motion-euclidean-for-the-warp-mode-in-ecc-image-alignment-method
    """
    print('\tAligning stacks')
    if save_as_uint16:
        s2_org = s2.copy()
        s2_corrected = np.zeros_like(s2)
        s1 = ff.rescale_to_uint8(s1) #Convert to unsigned integer 8-bit
        s2 = ff.rescale_to_uint8(s2) #Convert to unsigned integer 8-bit
    else:
        s2_corrected = np.zeros_like(s2)
    if (s1.dtype != np.uint8) or (s2.dtype != np.uint8):
        raise ValueError('s1 or s2 are not uint8!')
    
    #https://python.hotexamples.com/examples/cv2/-/findTransformECC/python-findtransformecc-function-examples.html
    # ff.pi(s1)
    # ff.pi(s2)
    im1 = np.median(s1, axis = 0).astype(np.float32)      
    im2 = np.median(s2, axis = 0).astype(np.float32)
   	# Find size of image1
    sz = im1.shape
  
    if not im1.shape == im2.shape:
            print("Input images have to be of the same size")
            return None    
        
   	# Define the motion model
    warp_mode = cv2.MOTION_TRANSLATION
   	
   	#Define the warp matrix
    warp_matrix = np.eye(2, 3, dtype=np.float32)
   
   	#Define the number of iterations
    number_of_iterations = 5000 #an integer between 5 and 5000
   
   	#Define correllation coefficient threshold
   	#Specify the threshold of the increment in the correlation coefficient between two iterations
    termination_eps = 1e-10#askfloat("Threshold", "Enter a number between 1e-10 and 1e-50",initialvalue=1e-10,minvalue=1e-50,maxvalue=1e-10)
   
   	#Define termination criteria
    criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations,  termination_eps)

   	#Run the ECC algorithm. The results are stored in warp_matrix.
    (cc, warp_matrix) = cv2.findTransformECC (im1,im2,warp_matrix, warp_mode, criteria)
 
    Z, Y, X = s2.shape
    for i in range(Z):
        if save_as_uint16:
            temporary_image = s2_org[i]

            # Use warpAffine for Translation, Euclidean and Affine
            aligned_image = cv2.warpAffine(temporary_image, warp_matrix, (sz[1], sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)

            s2_corrected[i] = aligned_image
        else:
            temporary_image = s2[i]

            # Use warpAffine for Translation, Euclidean and Affine
            aligned_image = cv2.warpAffine(temporary_image, warp_matrix, (sz[1], sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)
            
            s2_corrected[i] = aligned_image.astype(np.uint8)
    # print('3!!!')
    # ff.pi(s2_corrected)
    return s2_corrected


def match_size(I1, I2):
    """Matches the size of I1 and I2 in the x-dimension. Returns modified arrays"""
    if I1.shape[2] > I2.shape[2]: #First image is widest
        I1 = I1[:,:,0:-2]
    elif I1.shape[2] < I2.shape[2]: #Second image is widest
        I2 = I2[:,:,1:-1]
    if I1.shape != I2.shape:
        raise Exception('The dimensions of the 2 images are not the same after adjustment: '+str(I1.shape[1])+' and '+str(I2.shape[1]))
    return I1, I2

def split_image(I,A=0):
    """Returns two images that are half the image each in the x-direction
    I is in the dimension order of ZYX
    """
    if np.all(A == 0):
        w_I = I.shape[2]
        w_half_I = round(w_I/2)
        I_left = I[:,:,0:w_half_I]
        I_right = I[:,:,w_half_I+1:-1]
    else: #(x, y, width, height)    
        x1 = A[0][0]
        x2 = x1 + A[0][2]
        y1 = A[0][1]
        y2 = y1 + A[0][3]    
        x3 = A[1][0]
        x4 = x3 + A[1][2]
        y3 = A[1][1]
        y4 = y3 + A[1][3]   
        I_left = I[:,y1:y2,x1:x2]
        I_right = I[:,y3:y4,x3:x4]
    return I_left, I_right 

def show_RGB_based_on_two_images(I1,I2, title='Aligned', sz=10, rescale_individ=False, show=False, save=False, file_name=' ', frame_to_show = 0):
    if rescale_individ:
        #Rescales the pixel values of each image individually
        I1 = exposure.rescale_intensity(I1)
        I2 = exposure.rescale_intensity(I2)
    else:
        #Rescales the pixel values with shared min and max limits
        p_min = 5
        p_high = 99
        p_low_I1 = np.percentile(I1, p_min)
        p_high_I1 = np.percentile(I1, p_high)
        p_low_I2 = np.percentile(I2, p_min)
        p_high_I2 = np.percentile(I2, p_high)  
        p_low_min = np.min([p_low_I1,p_low_I2])  
        p_high_max = np.max([p_high_I1,p_high_I2])  
        I1 = exposure.rescale_intensity(I1, in_range=(p_low_min, p_high_max), out_range=np.uint8)
        I2 = exposure.rescale_intensity(I2, in_range=(p_low_min, p_high_max),out_range=np.uint8)
        im_RGB = np.transpose([I1, I2, np.zeros_like(I1)], (1, 2, 3, 0))

    if save: #save to file
        currentdir = os.path.dirname("__file__")
        dir = os.path.join(currentdir, 'polarisation\images', file_name)
        if not os.path.exists(dir):
            os.mkdir(dir)
        file_name = file_name+'_merged.tif'
        file_path = os.path.join(dir, file_name)
        print(im_RGB.shape)
        ti.write_tif_file(file_path, im_RGB, photometric='rgb')
    if show:
        #Adjust original image    
        fig = plt.figure(figsize=(sz, sz))  # create a figure object
        ax = fig.add_subplot(1, 1, 1)  # create an axes object in the figure
        ax.axis('off')
        ax.imshow(im_RGB[frame_to_show]) 
        plt.title(title)
        plt.show()  
    return im_RGB

    
#Is this function used?
def combine_into_stack(I1,I2):
    """
    Combines an split image into a stack. One of the images might be wider than the other and thus, this is adjusted for.
    """
    if I1.shape[1] > I2.shape[1]:
        #Left image is widest
        I1 = I1[:,0:-2]
    elif I1.shape[1] < I2.shape[1]:
        #Right image is widest
        I2 = I2[:,1:-1]
    if I1.shape != I2.shape:
        raise Exception('The dimensions of the 2 images are not the same after adjustment: '+str(I1.shape[1])+' and '+str(I2.shape[1]))
    
#     I_stacked = np.array([I1, I2])
    I_stacked = np.dstack((I1, I2))
    return I_stacked

#Is this function used?
def combine_into_stack_ZYX(I1,I2):
    """
    Combines an split image into a stack. One of the images might be wider than the other and thus, this is adjusted for.
    """
    if I1.shape[1] > I2.shape[1]:
        #Left image is widest
        I1 = I1[:,0:-2]
    elif I1.shape[1] < I2.shape[1]:
        #Right image is widest
        I2 = I2[:,1:-1]
    if I1.shape != I2.shape:
        raise Exception('The dimensions of the 2 images are not the same after adjustment: '+str(I1.shape[1])+' and '+str(I2.shape[1]))
    
#     I_stacked = np.array([I1, I2])
    I_stacked = np.array([I1, I2])
    return I_stacked


def load_made_polarization_files(v, settings_general, dir_exp, dic={}, frame_range = [-1,-1], only_I_tot=False, subtract_bg=True, rot90deg_after_split=True):
    if frame_range[0] == -1:
        frame_range = settings_general['frame_range']

    #Substring in the file name in order to find the right file    
    if 'prefix' in dic:
        prefix = dic['prefix']
    else:
        prefix = ''
    if 'sub_string' in dic:
        sub_string = dic['sub_string']
    else:
        sub_string = '0-0'
    if 'read_only_frame_subset' in dic:
        read_only_frame_subset = dic['read_only_frame_subset']
    else:
        read_only_frame_subset = False

    dir_pol = os.path.join(v.dir_path_file_folder, 'polarization_imgs')    
    if read_only_frame_subset:
        sub_string_perpen = 'perpen'+'_'+sub_string
        print(f'\tFinding files with prefix={prefix}, substring: {sub_string_perpen}')
        file_path_perpen = fp.find_file_name(dir_pol, sub_string=sub_string_perpen, prefix=prefix, file_ending = 'tif', return_path=True)
        file_path_parallel = file_path_perpen.replace("perpen", "parallel")
        file_path_tot = file_path_perpen.replace("perpen", "tot")    
    else:
        file_name_tot = f'tot_0-0_'+v.file_name0+'.tif'
        file_path_tot = os.path.join(dir_pol, file_name_tot)
        file_name_perpen = f'perpen_0-0_'+v.file_name0+'.tif'
        file_path_perpen = os.path.join(dir_pol, file_name_perpen)
        file_name_parallel = f'parallel_0-0_'+v.file_name0+'.tif'
        file_path_parallel = os.path.join(dir_pol, file_name_parallel)

    if not os.path.exists(file_path_tot):
        #Create the files (full stacks so will be slow)
        _, _, _ = load_and_split_imgs_optosplit(v.file_path, settings_general, dir_exp, crop=v.crop, crop2=v.crop2, save_to_file=True, frame_range=[0,0], subtract_bg = subtract_bg, rot90deg_after_split=rot90deg_after_split)
    
    print('xxx', frame_range)
    I_tot = ti.read_tif_file(file_path_tot, frame_range = frame_range)
    if not only_I_tot:
        I_perpen = ti.read_tif_file(file_path_perpen, frame_range = frame_range)
        I_parallel = ti.read_tif_file(file_path_parallel, frame_range = frame_range)
        return I_perpen, I_parallel, I_tot
    else:
        return I_tot

def save_polarization_files(v, frame_range, I_left, I_right_mod, I_tot, overwrite_saved_files):
    dir_save = os.path.join(v.dir_path_file_folder, 'polarization_imgs')
    if not os.path.exists(dir_save):
        os.mkdir(dir_save)
    #Save to two separate files:
    file_name1 = 'perpen_'+str(frame_range[0])+'-'+str(frame_range[-1])+'_'+v.file_name0+'.tif'
    file_path1 = os.path.join(dir_save, file_name1)
    file_name2 = 'parallel_'+str(frame_range[0])+'-'+str(frame_range[-1])+'_'+v.file_name0+'.tif'
    file_path2 = os.path.join(dir_save, file_name2)
    file_name_tot = 'tot_'+str(frame_range[0])+'-'+str(frame_range[-1])+'_'+v.file_name0+'.tif'
    file_path_tot = os.path.join(dir_save, file_name_tot)
    #Saves the files if they do not already exist or if they are to be overwritten.
    if not (os.path.exists(file_path1) & os.path.exists(file_path2)) or overwrite_saved_files:
        ti.write_tif_file(file_path1, I_left, photometric='minisblack')   
        ti.write_tif_file(file_path2, I_right_mod, photometric='minisblack') 
        ti.write_tif_file(file_path_tot, I_tot, photometric='minisblack')    
        os.startfile(dir_save) #Open file in explorer (only works for Windows)

def load_and_split_imgs_optosplit(file_path, settings_general, dir_exp, crop=[0, 0, 0, 0], crop2=[0, 0, 0, 0], save_to_file=False, frame_range=[0,0], overwrite_saved_files=False, subtract_bg = False, show_imgs_optosplit = False, align_manually=False, dx = 0, dy=0, lims_rescaling=[], rot90deg_after_split=False):
    """[Compares the alignment of the two sides frame-by-frame.]

    Args:
        file_path ([type]): [Contains a dual-split image (x-direction, 50/50 split)]
        crop (list, optional): [x, y, width, height]. Defaults to [0, 0, 0, 0]. 
        crop2 (list, optional): [x, y, width, height]. Defaults to [0, 0, 0, 0].
        save_to_file (bool, optional): [description]. Defaults to False.
        frame_range (list, optional): [description]. Defaults to [0,0].
        dir_save (str, optional): [description]. Defaults to ''.

    Raises:
        ValueError: [description]

    Returns:
        [I_left] (numpy array, uint8): [Left part of the image]
        [I_right_mod] (numpy array, uint8): [Right part of the image (modified so it is the same size as the left)]
    """
    if settings_general['show_imgs']:
        show_imgs_optosplit = True # for debugging

    save_as_uint16 = False
    if 'save_as_uint16' in settings_general:
        if settings_general['save_as_uint16']:
            save_as_uint16 = True

    tic = time.perf_counter()
    print("== Splitting the image in half and adjusting relative positions ==")

    #Read image
    file_path_copy = file_path.replace('\\', '.')
    file_path_copy = file_path_copy.replace('/', '.')
    # file_name0 = file_path_copy.split('.')[-2]

    #If the selected frames to read are all set to 0, read the full range in all dimensions. Otherwise, read only the selected range of dimensions.
    if file_path[-4:] == '.nd2':
        v = nd2_handling2.Video_waves(file_path, read_img_directly=True, frame_range = frame_range, dir_exp = dir_exp, subtract_bg=subtract_bg)  
        if hasattr(v,'file_path_bg'):
            file_path_bg = v.file_path_bg
        img = v.img
    if file_path[-5:] == '.tiff' or file_path[-4:] == '.tif':
        raise ValueError('file type must be nd2...')

    if subtract_bg & hasattr(v,'file_path_bg'):
        img = ff.subtract_uneven_bg_from_illumination(img, file_path_bg, ignore_black_pixels=True)

    if save_as_uint16:
        img = img.astype('uint16')
    else:
        #If the light intensity is low, re-scale the image unevenly to uint8, saving the pixel value resolution
        # I_max = np.max(img)
        I_max = np.percentile(img,99.9999)
        # ff.show_img(img[0], title_text = 'Before split', cmap='gray', sz=5, fsz = 15)
        I_threshold = 20000 #Threshold for when to do re-scaling to uint8 with lower upper I limit
        print(f'\tMaximum intensity of the image stack: {I_max:.0f}, threshold: {I_threshold}')        
        # print('\nBefore re-scaling to uint8')
        # ff.pi(img)
        if I_max < I_threshold:
            img = ff.rescale_to_uint8(img, lims=[0,I_max]) #Convert to unsigned integer 8-bit  
        else:
            img = ff.rescale_to_uint8(img) #Convert to unsigned integer 8-bit  
        # else:    
        #     img = ff.rescale_to_uint8(img) #Convert to unsigned integer 8-bit  
        # print('\nAfter re-scaling to uint8')
        # ff.pi(img)             
    #Raise an exception if the image has the wrong dimensions.
    if len(img.shape) != 3:
        raise ValueError(f"The image must be 3D and not {len(img.shape)}D")

    #Split image into two halves (3D)
    I_left, I_right = split_image_in_half(img)
    if show_imgs_optosplit:
        ff.show_img(I_left[0], title_text = 'After split (Left image)', cmap='gray', sz=5, fsz = 15)
        ff.show_img(I_right[0], title_text = 'After split (Right image)', cmap='gray', sz=5, fsz = 15)
        
        print('shape left image: ', I_left.shape)
        print('shape right image: ', I_right.shape)

    # #Crop images 1 (to get rid of the black strip in between and next to the frames)
    if not all(i == 0 for i in crop):        
        I_right = ff.crop_image(I_right, crop)
        I_left = ff.crop_image(I_left, crop)

        if show_imgs_optosplit:
            ff.show_img(I_left[0], title_text = 'After 1st cropping\n(Left image)', cmap='gray', sz=5, fsz = 15)
            ff.show_img(I_right[0], title_text = 'After 1st cropping\n(Right image)', cmap='gray', sz=5, fsz = 15) 

    #Align images with my script
    if align_manually:
        I_right_mod = align_two_stacks_manually(I_left, I_right, save_as_uint16, dx=dx, dy=dy) #Right image aligned to left
    else:
        I_right_mod = align_two_stacks(I_left, I_right, save_as_uint16) #Right image aligned to left
    if show_imgs_optosplit:
        print('shape left image: ', I_left.shape)
        print('shape right image: ', I_right.shape)
        print('shape right image after mod: ', I_right_mod.shape)

    if show_imgs_optosplit:
        ff.show_img(I_right_mod[0], title_text = 'After alignment\n(Right image)', cmap='gray', sz=5, fsz = 15)

    # #Crop images 2 (To get rid of the black border in case the right side moved)
    if not all(i == 0 for i in crop2):
        print(f'\tCropping image according to: {crop2}')
        I_right_mod = ff.crop_image(I_right_mod, crop2)
        I_left = ff.crop_image(I_left, crop2)
        print('shape left image: ', I_left.shape)
        print('shape right image: ', I_right.shape)
        print('shape right image after mod: ', I_right_mod.shape)
        if show_imgs_optosplit:
            ff.show_img(I_left[0], title_text = 'After 2nd cropping (Left image)', cmap='gray', sz=5, fsz = 15)              

            ff.show_img(I_right_mod[0], title_text = 'Right image after alignment and 2nd cropping', cmap='gray', sz=5, fsz = 15) 
            dst = ff.blend_images(I_left[0], I_right_mod[0])
            ff.show_img(dst, sz=20, cmap='gray', title_text='Overlay after cropping, right mod on left')

    if rot90deg_after_split and (not save_as_uint16):
        #Rotate images
        print('\tRotating both left and right side 90 degrees.')
        I_left = ff.transform_img(I_left, angle=90, transform_mode='rot')
        I_right_mod = ff.transform_img(I_right_mod, angle=90, transform_mode='rot')

    if show_imgs_optosplit:
        im_merged = show_RGB_based_on_two_images(I_left,I_right_mod)
        ff.show_img(im_merged[0], title_text = 'Red and Green', cmap='gray', sz=5, fsz = 15)
        ff.show_img(im_merged[0,130:190,130:190], title_text = 'Red and Green (zoomed-in)', cmap='gray', sz=5, fsz = 15)

    #Create I_tot
    I_tot = np.mean( np.array([I_left, I_right_mod]), axis=0)

    if not save_as_uint16:
        I_tot = ff.rescale_to_uint8(I_tot) #Convert to unsigned integer 8-bit  

    if save_to_file:
        save_polarization_files(v, frame_range, I_left, I_right_mod, I_tot, overwrite_saved_files)
    toc = time.perf_counter()
    print(f"\tSplit and adjusted the two-sided image in {toc - tic:0.2f} seconds")
    return I_left, I_right_mod, I_tot

def split_image_in_half(I):
    """Splits a numpy array into two halves"""
    w_I = I.shape[2]
    w_half_I = round(w_I/2)
    I_left = I[:,:,0:w_half_I]
    I_right = I[:,:,w_half_I+1:-1]
    I_left, I_right = match_size(I_left, I_right)
    return I_left, I_right

def calcAnisotropy(I1_img, I2_img, range=(0,1/3), p=[-0.1, 0.1]):
    """Calculates and return the Anisotropy as:
        I_aniso = (I1-I2) / (I1 + 2*I2)
        I1 is the image parallel to the flow direction, usually the right side of the image
        I2 is the image perpendicular to the flow direction, usually the left side of the image
    """

    tic = time.perf_counter()
    print('Calculating the Anisotropy')
    I_aniso = calc_anisotropy_raw(I1_img, I2_img)
    I_aniso = exposure.rescale_intensity(I_aniso, in_range=(p[0], p[1]), out_range=range).astype(np.float16)
    toc = time.perf_counter()
    print(f"\tCalculated the anisotropy in {toc - tic:0.2f} seconds")
    return I_aniso

def calc_p_ratio(I_parallel, I_perpen):
    "Calculates the emission polarization ratio, P. See the paper van Mameren, Vermeulen et al. 2018 for more information."
    I_parallel = I_parallel.astype(float)
    I_perpen = I_perpen.astype(float)    
    P = (I_parallel - I_perpen) / (I_parallel+I_perpen)
    return P

def calc_anisotropy_raw(I1_img, I2_img):
    """Calculates and return the raw Anisotropy as:
        I_aniso = (I1-I2) / (I1 + 2*I2)
        I1 is the image parallel to the flow direction, usually the right side
        I2 is the image perpendicular to the flow direction, usually the left side
    """
    raise ValueError("Polarization calculation type: 'anisotropy' is currently not eligible to minimize error")
    I1 = I1_img.astype(float)
    I2 = I2_img.astype(float)
    I_aniso = (I1-I2) / (I1 + 2*I2)   
    I_aniso = np.nan_to_num(I_aniso, copy=True, nan=0.0, posinf=None, neginf=None) 
    return I_aniso

from matplotlib.colors import hsv_to_rgb
def convert_to_hsv_space(I_tot, h, rescale=True, in_intensity_range=(0,0), ):
    """
    Convert to a HSV-image, where the hue represents the polarisation fluorescence anisotropy, r,  
    or polarization emission ratio, P, ranging from 0 (red) to 0.33 (green). 

    The HSV image is then converted to a RGB image
    
    Info about HSV
    HSV-image values in matplotlib range from 0-1 in all three channels: hue, saturation and value. 
    - Hue values range from 0 to 360. We want values from red (0) to green (120)
    - Value (Brightness), should be a value between 0 and 1
    - Saturation. We don't use saturation to display information.
    
    """
    tic = time.perf_counter()
    print('Converting to HSV space')
  
    #Saturation, keep it at maximum (=1)
    s = np.ones_like(h, dtype=np.uint8)
    
    #Adjust the min and max value for the LUT to enhance the contrast:    
    if rescale == True:
        v = exposure.rescale_intensity(I_tot, in_range=in_intensity_range, out_range=(0, 1)).astype(np.float16) #value 
    else:
        v = I_tot.astype(np.float16) 

    print(f'Value part of HSV (I_tot) based on following brightness/constrast limits: {in_intensity_range}')
    #HSV image, transpose it to convert the shape to ZXYC
    if I_tot.ndim == 3:
        im_hsv = np.transpose([h, s, v], (1, 2, 3, 0))
    elif I_tot.ndim == 2:
        im_hsv = np.transpose([h, s, v], (1, 2, 0))

    del I_tot #clear memory
    del h, s, v #clear memory
    toc = time.perf_counter()
    print(f'\tconverted to hsv after {toc - tic:0.2f} seconds')

    #Convert from HSV to RGB
    im_rgb = hsv_to_rgb(im_hsv)
    del im_hsv
    toc = time.perf_counter()
    print(f'\tconverted to rgb after {toc - tic:0.2f} seconds')

    #Rescale the RGB image to values between 0 and 255
    im_rgb = exposure.rescale_intensity(im_rgb,out_range=(0,255)).astype(np.uint8)
    print('\trescaled and converted to uint8')

    toc = time.perf_counter()
    print(f"\tConverted to HSV space in {toc - tic:0.2f} seconds")
    return im_rgb

import mpl_toolkits.axes_grid1.anchored_artists as mpl_aa
import mpl_toolkits.axes_grid1.inset_locator as mpl_il
import matplotlib.colors
def add_colourwheel(axis, loc='upper left', size=.12, **wheel_kw):
    """ 
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger    
    Add a colourwheel to a plot. 
    Use `wheel_kw` to choose alpha, saturation, etc of the wheel image
    """

    wheel_im = _generate_wheel(**wheel_kw)

    inset = mpl_il.inset_axes(
        axis, loc='upper left', width=size, height=size
    )
    inset.axis('off')
    inset.imshow(wheel_im)

def add_colourwheel_in_place(im, pad, resolution):
    """ 
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger
    Add a colourwheel to an np.array
    Use `wheel_kw` to choose alpha, saturation, etc of the wheel image
    """
    wheel_im = _generate_wheel(resolution=resolution)
    Y, X, _ = wheel_im.shape
    wheel_rgb = wheel_im[..., 0:3]
    wheel_alpha = wheel_im[..., 3].astype(int).reshape(Y, X, 1)

    im = im.copy()
    im[pad:pad+Y, pad:pad+X] = wheel_alpha * wheel_rgb + \
        (1-wheel_alpha) * im[pad:pad+Y, pad:pad+X]
    return im


def _generate_wheel(alpha=1, sat=1, val=.9, resolution=100):
    """
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger
    """
    x = np.linspace(-1, 1, resolution)
    y = np.linspace(-1, 1, resolution)

    xx, yy = np.meshgrid(x, y)
    im_hsv = np.zeros((*xx.shape, 3))

    rr = np.sqrt(xx**2 + yy**2)  # Distance from origin
    tt = np.arctan2(yy, -xx)     # Angle of pixel vector to x axis

    disc = rr < 1  # Image that is 0 if pixel outside unit circle, 1 otherwise

    im_hsv[..., 0] = 1 - (disc*tt / np.pi + 1) % 1  # Hue
    im_hsv[..., 1] = disc*rr * sat                # Saturation
    im_hsv[..., 2] = disc * val                   # Brightness

    # Convert to RGB
    im_hsv = matplotlib.colors.hsv_to_rgb(im_hsv)
    im_rgb = np.zeros((*xx.shape, 4))
    im_rgb[..., 0:3] = im_hsv

    # Add alpha channel
    # alpha should be 0 outside the disk, otherwise the image will be black
    im_rgb[..., 3] = disc * alpha
    return im_rgb

def add_circle_sector(axis, loc='upper left', size=.12, **wheel_kw):
    """ 
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger    
    Add a colourwheel to a plot. 
    Use `wheel_kw` to choose alpha, saturation, etc of the wheel image
    """

    wheel_im = _generate_circle_sector(**wheel_kw)

    inset = mpl_il.inset_axes(
        axis, loc='upper left', width=size, height=size
    )
    inset.axis('off')
    inset.imshow(wheel_im)

def add_circle_sector_in_place(im, pad, resolution):
    """ 
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger
    Add a colourwheel to an np.array
    Use `wheel_kw` to choose alpha, saturation, etc of the wheel image
    """
    wheel_im = _generate_circle_sector(resolution=resolution)
    Y, X, _ = wheel_im.shape
    wheel_rgb = wheel_im[..., 0:3]
    wheel_alpha = wheel_im[..., 3].astype(int).reshape(Y, X, 1)
    wheel_rgb = exposure.rescale_intensity(wheel_rgb, out_range=(0,255)).astype(np.uint8)
    im = im.copy()
    im[pad:pad+Y, pad:pad+X] = wheel_alpha * wheel_rgb + \
        (1-wheel_alpha) * im[pad:pad+Y, pad:pad+X]
    return im

def _generate_circle_sector(alpha=1, sat=1, val=1, resolution=100):
    """
    Copied from https://github.com/wduverger/tegen/blob/main/tegen/polarisation.py
    Author: Wouter Duverger
    """
    x = np.linspace(0, 1, resolution) #changed range from 0 to 1
    y = np.linspace(0, 1, resolution) #changed range from 0 to 1

    xx, yy = np.meshgrid(x, y)
    im_hsv = np.zeros((*xx.shape, 3))

    rr = np.sqrt(xx**2 + yy**2)  # Distance from origin
    tt = np.arctan2(yy, -xx)     # Angle of pixel vector to x axis

    disc = rr < 1  # Image that is 0 if pixel outside unit circle, 1 otherwise

    #im_hsv[..., 0] = 1 - (disc*tt / np.pi + 1) % 1  # Hue
    im_hsv[..., 0] = (1 - (disc*tt / np.pi))*240/360 # Hue (updated)
    im_hsv[..., 1] = disc*rr * sat                # Saturation
    im_hsv[..., 2] = disc * val                   # Brightness

    # Convert to RGB
    im_hsv = matplotlib.colors.hsv_to_rgb(im_hsv)
    im_rgb = np.zeros((*xx.shape, 4))
    im_rgb[..., 0:3] = im_hsv

    # Add alpha channel
    # alpha should be 0 outside the disk, otherwise the image will be black
    im_rgb[..., 3] = disc * alpha
    return im_rgb

import pandas as pd
import os
        
def show_subplots_anisotropy(frame, I_parallel, I_perpen, I_tot, im_merged, im_rgb, I_aniso, v, aniso_hsv_cbar_range=[-0.1, 0.1], p_range_hist=[-0.2,0.2], sz=[5,2], fsz=15, sub_folder='', pol_type_text = 'P'):
    img_hsv = copy.deepcopy(im_rgb[frame])
    

    # if sub_folder == 'pixelated':
    #     add_aniso_cbar = -1
    # else:
    #     add_aniso_cbar = 3
    #     img_hsv = add_circle_sector_in_place(img_hsv, pad=10, resolution=50)

    from skimage import exposure
    if pol_type_text == 'r':
        str_pol = '$r=(I_{\parallel} - I_{\perp} ) / (I_{\parallel} + 2\cdot I_{\perp})$'
    elif pol_type_text == 'P':
        str_pol = r'$P=(I_{\parallel} - I_{\perp} ) / (I_{\parallel} + I_{\perp})$'
    else:
        raise ValueError('Wrong pol_type_text input')
    I_list = [I_parallel[frame], 
              I_perpen[frame], 
              exposure.rescale_intensity(I_tot[frame], out_range=(0,255)).astype(np.uint8), 
              im_merged[frame],
              I_aniso[frame],
              img_hsv              
             ]
    titles = ['$I_{\perp}$',
              '$I_{\parallel}$',
              '$I_{\perp} + I_{\parallel}$', 
              '$I_{\perp}$ red\n$I_{\parallel}$ green',
              str_pol, 
              f'{pol_type_text} gives color\n and I gives brightness\n {pol_type_text} range: [{aniso_hsv_cbar_range[0]}, {aniso_hsv_cbar_range[1]}]']
    ff.plot_sub_plots(I_list, titles, n_cols=3, sz=sz, fsz = fsz, cmap='gray', suptit=f'{v.p} mbar, frame #{frame}')
  
    #Histogram of the pixel intensities of the two images
    fig, ax = plt.subplots(figsize=(15,4))
    plt.style.use('ppt') 
    I_p = I_perpen[frame].ravel()
    I_p = I_p[I_p != 0] 
    I_par = I_parallel[frame].ravel()
    I_par = I_par[I_par != 0] 
    plt.hist(I_p,bins=256, color='r', range=(0,max(I_perpen.ravel())), label='$I_{\perp}$', alpha=0.5)
    plt.hist(I_par,bins=256, color='g', range=(0,max(I_perpen.ravel())), label='$I_{\parallel}$', alpha=0.5)
    plt.title(f'Histogram of pixel intensities, {v.p} mbar, frame #{frame}', fontsize=fsz)
    plt.legend(loc='upper right')
    plt.xlabel('pixel values')
    plt.ylabel('count')  
    # plt.style.use('general') 
    plt.show()

    print(f'Information about Image_{pol_type_text}[{frame}]')
    ff.pi(I_aniso[frame])  

    A = I_aniso[frame].ravel()
    # A = A[A != 0]  
    #Anisotropy Histogram
    fig, ax = plt.subplots(figsize=(15,4))
    plt.style.use('ppt')
    plt.title(f'Histogram of {pol_type_text}'+r' $=(I_{\parallel} - I_{\perp} ) / (I_{\parallel} + 2\cdot I_{\perp})$'+f' {v.p} mbar, frame #{frame}')
    y, x, _ = plt.hist(A,bins=256, color='r', range=p_range_hist)
    plt.xlabel(pol_type_text)
    plt.ylabel('count')   
    plt.ylim([0,np.max(y)/3])      
    p = ff.calc_percentiles(I_aniso[frame], min=2, max=98)
    plt.axvline(p[0], color='b', label='2% percentile')
    plt.axvline(p[1], color='b', label='98% percentile')
    plt.legend(loc='upper right')
    plt.show()

def plot_quiver_vector_field(img1, img2, ax, sz=5, autoscaling_on=True, p_lim = [2,98], n_vectors_per_row=10, cmap='inferno'):
    """Plots a vector field (arrows) on top of an image corresponding to
    the radius of the two vectors corresponding to the pixel value of each image.
    Autoscales by default.

    Args:
        img1 ([type]): [description]
        img2 ([type]): [description]
        frame ([type]): [description]
        sz (int, optional): [description]. Defaults to 5.
        autoscaling_on (bool, optional): [description]. Defaults to True.
        p_lim (list, optional): [description]. Defaults to [2,98].
        n_vectors_per_row (int, optional): [description]. Defaults to 10.
        cmap (str, optional): [description]. Defaults to 'viridis'.
    """
    r = np.sqrt(img1**2 + img2**2) #radius
    w_img=img1.shape[-1] #Image width
    #Calculate the percentiles of the pixel values for brightness/contrast adjustment
    p = ff.calc_percentiles([img1, img2], min=p_lim[0], max=p_lim[1])
    # print(f'Percentiles ({p_lim[0]}, {p_lim[1]}) equals pixel values ({p[0]},{p[1]})')
    if autoscaling_on:
        # scaling_factor = round(w_img/n_vectors_per_row)*p_lim[1]/100
        scaling_factor = round(w_img/n_vectors_per_row)*0.75
        # print(f'Autoscales with factor = {scaling_factor}')  

    x,y = np.meshgrid(range(0,img1.shape[1]), range(0,img1.shape[0]))
    # u = img1[frame]*-1 #vector field in one direction
    # v = img2[frame] #vector field in the other direction
    u = ((img1 - p[0])/(p[1] - p[0]))*scaling_factor
    u=u*-1
    v = ((img2 - p[0])/(p[1] - p[0]))*scaling_factor        

    # fig, ax = plt.subplots(figsize=(sz,sz))
    # Every j:rd point in each direction.
    j=round(w_img/n_vectors_per_row)
    skip = (slice(None, None, j), slice(None, None, j))
    ax.quiver(x[skip],y[skip],u[skip],v[skip], angles='xy', scale_units='xy', scale=1, color='white') #, width=0.005*0.12
    #ax.set_ylim(ax.get_ylim()[1], ax.get_ylim()[0])
    # im = ax.imshow(np.flipud(r[frame]), extent=[x.min(), x.max(), y.min(), y.max()], alpha=0.5)
    im = ax.imshow(r, alpha=1, cmap=cmap)
    #ax.set(aspect=1, title=f'frame #{frame}')
    # plt.axis('off')

from matplotlib import patches
def draw_ellipse(x0,y0,w,h,a, ax):
    e = patches.Ellipse((x0, y0), w, h,
                         angle=a, linewidth=1, fill=False, zorder=2)
    ax.add_patch(e)

def draw_ellipse_field(img1,img2,img_show, ax, n_vectors_per_row=5, scaling_factor=10, sz=5,title='', p_lim = [2,98],autoscaling_on = True):    
    
    w_img=img1.shape[-1]
    h_img=img1.shape[-2]
    n_x = n_vectors_per_row
    n_y = n_vectors_per_row
    angle=0
    p = ff.calc_percentiles([img1, img2], min=p_lim[0], max=p_lim[1])
    # print(f'Percentiles ({p_lim[0]}, {p_lim[1]}) equals pixel values ({p[0]},{p[1]})')
    if autoscaling_on:
        scaling_factor = round(w_img/n_x)*p_lim[1]/100
        # print(f'Autoscales with factor = {scaling_factor}')

    # fig, ax = plt.subplots(figsize=(sz,sz))
    x = np.linspace(0,w_img,n_x)
    y = np.linspace(0,h_img,n_y)
    for i in range(x.shape[0]-1):
        for j in range(y.shape[0]-1):
            x0=x[i]+x[1]/2
            y0=y[j]+y[1]/2
            a = ((img1[round(y0),round(x0)]-p[0])/(p[1]-p[0]))*scaling_factor
            b = ((img2[round(y0),round(x0)]-p[0])/(p[1]-p[0]))*scaling_factor              
            #print(f'x={round(x0)}, \ty={round(y0)}, \t a={np.round(a,1)}, \t b={np.round(b,1)}, \tangle={angle}')
            draw_ellipse(x0,y0,a,b,angle, ax)
    im = ax.imshow(img_show, alpha=1, cmap='inferno')

def show_histogram_I_tot(I_tot, in_intensity_range):
    """Shows the histogram of the Total Intensity For all pixels in all frames in the image stack I_tot
    Also displays the intensity limits (in_intensity_range) as vertical dashed lines

    Args:
        I_tot (_type_): _description_
        in_intensity_range (_type_): _description_
    """
    fig, ax = plt.subplots(figsize=(15,4))
    plt.style.use('ppt') 
    I = I_tot.ravel()
    I = I[I != 0] 
    n_black = len(I_tot) - len(I)
    plt.hist(I,bins=256, color='b', range=(0,max(I_tot.ravel())), label='$I_{\perp}$')
    if n_black > 0:
        plt.title('Histogram of I_tot pixel intensities\n'+f'(n black pizels = {n_black})')
    else:
        plt.title(f'Histogram of I_tot pixel intensities')
    plt.xlabel('pixel values')
    plt.ylabel('count')  
    plt.axvline(in_intensity_range[0], linestyle='--', color='k', label=f'I1 = {in_intensity_range[0]}')
    plt.axvline(in_intensity_range[1], linestyle='--', color='k', label=f'I2 = {in_intensity_range[1]}')
    plt.legend(loc='upper right')
    plt.show()

def calc_aniso_and_show_output(v, settings_general, d_vid, dic, frame=0, sub_folder='', p_range=(-0.2,0.2), range_hsv_out=(0,1/3), I_tot_percentiles=(0.5,99.5), I_tot_limits=(0,0), average_over_frames=False):
    """[summary]
    Either set the limits of the LUT for the total intensity using percentiles
    with 'I_tot_percentiles' or manually set the limits with 'I_tot_limits'


    Args:
        v ([type]): [description]
        settings_general ([type]): [description]
        frame (int, optional): [description]. Defaults to 0.
        sub_folder (str, optional): [description]. Defaults to ''.
        p : Range for the anisotropy to convert to color
    """
    #Substring in the file name in order to find the right file    
    if 'prefix' in dic:
        prefix = dic['prefix']
    else:
        prefix = 'dz_pix_compressed'
    if 'sub_string' in dic:
        sub_string = dic['sub_string']
    else:
        sub_string = '0-0'
    if 'rot90deg_after_split' in dic:
        rot90deg_after_split = dic['rot90deg_after_split']
    else:
        rot90deg_after_split = False
    if 'read_only_frame_subset' in dic:
        read_only_frame_subset = dic['read_only_frame_subset']
    else:
        read_only_frame_subset = False

    if read_only_frame_subset:
        extra_file_name_save_text = sub_string
    # mode = settings_general['mode']

    polarization_calc_type = settings_general['polarization_calc_type']

    # raise ValueError('Update this function.')
    print(f'  - Reading files, calculating {polarization_calc_type} \n\tv: File #{v.file_nbr}, frame range: {v.frame_range}, p = {v.p} mbar, frame rate = {v.frame_rate} fps, cropping area = {v.crop}')
    if sub_folder == 'pixelated': #pixelated image
        #Find file paths
        sub_string_perpen = 'perpen'+'_'+sub_string
        print(f'\tFinding files with prefix={prefix}, substring: {sub_string_perpen}')
        file_path_perpen = fp.find_file_name(v.dir_pixelated_files, sub_string=sub_string_perpen, prefix=prefix, file_ending = 'tif', return_path=True)
        file_path_parallel = file_path_perpen.replace("perpen", "parallel")
        file_path_tot = file_path_perpen.replace("perpen", "tot")

        #Read files
        I_parallel = ti.read_tif_file(file_path_parallel, frame_range = v.frame_range)
        I_perpen = ti.read_tif_file(file_path_perpen, frame_range = v.frame_range)
        I_tot = ti.read_tif_file(file_path_tot, frame_range = v.frame_range)
    else:#Non-pixelated, raw image
        I_perpen, I_parallel, I_tot = load_made_polarization_files(v, settings_general, d_vid['dir_exp'], dic=dic, frame_range = v.frame_range, rot90deg_after_split=rot90deg_after_split)

    if average_over_frames:
        I_tot=np.mean(I_tot, axis=(0))
        I_parallel=np.mean(I_parallel, axis=(0))
        I_perpen=np.mean(I_perpen, axis=(0))

    # If the limits have been set in I_tot_lims, use those. 
    # If not, calculate limits based on percentages of the pixel values
    if I_tot_limits[1] == 0:        
        in_intensity_range = ff.calc_percentiles(I_tot, min=I_tot_percentiles[0], max=I_tot_percentiles[1])
    else:
        in_intensity_range = I_tot_limits

    if settings_general['show_imgs']:
        show_histogram_I_tot(I_tot, in_intensity_range)

    #Calculate the fluorescence anisotropy or the polarization emission ratio
    if polarization_calc_type == 'p_ratio':
        I_pol_raw = calc_p_ratio(I_parallel, I_perpen)
        pol_type = 'p_ratio_raw'
        pol_type_text = 'P'
    else:
        raise ValueError('Wrong input of polarization_calc_type in settings_general')
    ff.pi(I_pol_raw)
    print(f'Range in {pol_type} for rescaling to color: {p_range}')
    I_pol_rescaled = exposure.rescale_intensity(I_pol_raw, in_range=(p_range[0], p_range[1]), out_range=range_hsv_out).astype(float)
    print(f'shape of {pol_type}:  {I_pol_raw.shape}')

    #Convert to the HSV-format and then into RGB-based arrays for saving into image files
    im_rgb = convert_to_hsv_space(I_tot, I_pol_rescaled, in_intensity_range=in_intensity_range)
    
    if settings_general['show_imgs'] and len(I_perpen.shape) == 3:
        im_merged = show_RGB_based_on_two_images(I_perpen,I_parallel)

        #Show an array of images (both original and modified)
        show_subplots_anisotropy(frame, I_parallel, I_perpen, I_tot, im_merged, im_rgb, I_pol_raw, v=v, sz=[15,9], aniso_hsv_cbar_range=p_range, sub_folder=sub_folder, pol_type_text=pol_type_text, p_range_hist=p_range)

    if settings_general['save_imgs']:
        #Save HSV image to file
        save_hsv_img(settings_general, v, im_rgb,p_range, extra_file_name_save_text=extra_file_name_save_text)

        #Save anistropy/p-ratio to file
        save_pol_raw_img(settings_general, v, I_pol_raw, p_range, pol_type=pol_type)

def save_hsv_img(settings_general, v, im_rgb,p, extra_file_name_save_text=''):
    """[Save HSV image to file, split into multiple file if requested]

    Args:
        settings_general ([type]): [description]
        v ([type]): [description]
        im_rgb ([type]): [description]
        p ([type]): [description]
    """

    file_name_hsv = str(v.frame_range[0])+'-'+str(v.frame_range[-1])+'_'+v.file_name0 +'_hsv_p'+str(p[1])+'.tif'
    file_name_hsv = file_name_hsv.replace('perpen_','')
    if len(extra_file_name_save_text) > 0:
        file_name_hsv = extra_file_name_save_text + '_' + file_name_hsv
    if 'hsv_split_file' in settings_general:
        if settings_general['hsv_split_file']:
            file_name_hsv = 'split_'+file_name_hsv
    if settings_general['mode_pixelation'] == 'dead_zone':
        file_name_hsv = 'dz_pix_'+file_name_hsv

    file_path_hsv = os.path.join(v.dir_hsv, file_name_hsv) 
    ti.write_tif_file(file_path_hsv, im_rgb, photometric='rgb')
    if settings_general['open_saved_file_dir_after_script']:
        os.startfile(v.dir_hsv) #Open file in explorer (only works for Windows)

def save_pol_raw_img(settings_general, v, im_rgb,p, pol_type = '_p_ratio_raw'):
    file_name_hsv = str(v.frame_range[0])+'-'+str(v.frame_range[-1])+'_'+v.file_name0 +pol_type+'.tif'
    file_name_hsv = file_name_hsv.replace('perpen_','')
    if 'hsv_split_file' in settings_general:
        if settings_general['hsv_split_file']:
            file_name_hsv = 'split_'+file_name_hsv
    if settings_general['mode_pixelation'] == 'dead_zone':
        file_name_hsv = 'dz_pix_'+file_name_hsv

    file_path_hsv = os.path.join(v.dir_hsv, file_name_hsv) 
    ti.write_tif_file(file_path_hsv, im_rgb, photometric='minisblack')
    if settings_general['open_saved_file_dir_after_script']:
        os.startfile(v.dir_hsv) #Open file in explorer (only works for Windows)

def split_files_calc_aniso_and_show_output(v,settings_general, d_vid, n_frames_limit, sub_folder, dic, p_range = (-0.2,0.2), I_tot_percentiles=(0.5,99.5), I_tot_limits=(0,0)):
    n_subimgs = math.ceil(settings_general['frame_range'][-1]/n_frames_limit)+1            
    settings_general['hsv_split_file'] = True
    ranges = np.round(np.linspace(settings_general['frame_range'][0],settings_general['frame_range'][-1],n_subimgs))
    ranges = np.array(ranges,dtype=int)
    print(f'\n\t Splitting up the processing into {len(ranges)} substacks of {ranges[1]-ranges[0]+1} frames to reduce memory load:')
    for i in range(0,len(ranges)-1): print(f'\t\t{ranges[i]}-{ranges[i+1]}')
    for k in range(len(ranges)-1):
        frame_start = ranges[k]
        frame_end = ranges[k+1] 
        if k == len(ranges):
            frame_end = frame_end + 2
            print('LAST! ')
        settings_general['frame_range'] = range(frame_start,frame_end)
        v.frame_range = range(frame_start,frame_end)
        print('frame_start',frame_start,'frame_end',frame_end)

        print('not showing result')
        calc_aniso_and_show_output(v, settings_general, d_vid, dic, frame=v.frame_to_display, sub_folder=sub_folder, p_range = p_range,I_tot_percentiles=I_tot_percentiles, I_tot_limits=I_tot_limits)

def split_imgs_polarization_mode(file_path, settings_general, d):
    """Split the optosplit (dual-polarization) images into separate files (left, right and mean of both left and right)"""

    if 'd_vid' in d:
        d_vid = d['d_vid']        
    else:
        raise ValueError('Could not find d_vid in the dictionary')
    if 'df_vid' in d:
        df_vid = d['df_vid']        
    else:
        raise ValueError('Could not find df_vid in the dictionary')

    #For some videos, e.g. some 100x videos, I forgot to get a BG illumination heterogeneity image. For these images, I skip adding the file path to the background image.    
    if 'subtract_bg' in d:
        subtract_bg = d['subtract_bg'] 
    else:    
        subtract_bg = True

    if 'subtract_bg' in d_vid:
        if d_vid['subtract_bg'] == 'False' or d_vid['subtract_bg'] == False:
            subtract_bg = False
            print('Skipping background subtraction as directed from the experiment spreadsheet')

    if 'align_manually' in d:
        align_manually = d['align_manually'] 
    else:    
        align_manually = False
    if 'dy' in d:
        dy = d['dy'] 
    else:    
        dy = 0
    if 'dx' in d:
        dx = d['dx'] 
    else:    
        dx = 0
    if 'rot90deg_after_split' in d:
        rot90deg_after_split = d['rot90deg_after_split']
    else:
        rot90deg_after_split = False

    #Cropping of both side frames before alignment of the two sides
    if 'crop' in d:
        crop = d['crop']
    else:
        crop = [0, 0, 0, 0]

    #Cropping of eachafter alignment of the two sides
    if 'crop2' in d:
        crop2 = d['crop2']
    else:
        crop2 = [0, 0, 0, 0]
    dir_exp = d_vid['dir_exp']

    if 'load_crop_values_from_spreadsheet' in d:
        if d['load_crop_values_from_spreadsheet']:            
            crop = ap.extract_crop_values_From_df(df_vid) #Crop settings of the raw image
            crop2 = ap.extract_crop_values_From_df(df_vid, str='crop_area2')
            print(f'\tLoaded crop values from spreadsheet: crop: {crop}, crop2: {crop2}')
    
    _, _, _ = load_and_split_imgs_optosplit(file_path, settings_general, dir_exp, crop=crop, crop2=crop2, save_to_file=settings_general['save_imgs'], frame_range=settings_general['frame_range'], subtract_bg = subtract_bg, align_manually=align_manually, dx=dx, dy=dy,rot90deg_after_split=rot90deg_after_split)
    
    return {}

import nd2_handling2
def hsv_image_generation_v2(file_path, settings_general_org, d):
    settings_general=settings_general_org.copy() #copy because settings_general is modified
    #Maximum number of frames to process in a single session (to minmimize too high CPU load at the same time)
    if 'n_frames_limit' in d:
        n_frames_limit = d['n_frames_limit']
    else:
        n_frames_limit = 500
    if 'df_vid' in d:
        df_vid = d['df_vid']
    else:
        raise ValueError('Could not find df_vid in the dictionary')
    if 'd_vid' in d:
        d_vid = d['d_vid']        
    else:
        raise ValueError('Could not find d_vid in the dictionary')

    #Limits for polarization emission LUTs (which will turn into hue in the HSV-model)
    if 'p_range' in d:
        p_range = d['p_range']        
    else:
        p_range = (-0.2,0.2)

    #Limits for total intensity LUTs (which till turn into value in the HSV-model)
    if 'I_tot_percentiles' in d:
        I_tot_percentiles = d['I_tot_percentiles']
    else:
        I_tot_percentiles=(0.5,99.5)

    #Limits for total intensity LUTs (which till turn into value in the HSV-model)
    if 'I_tot_limits' in d:
        I_tot_limits = d['I_tot_limits']
    else:
        I_tot_limits=(0,0)

    if 'average_over_frames' in d:
        average_over_frames = d['average_over_frames']
    else:
        average_over_frames=False
    if 'read_only_frame_subset' in d:
        read_only_frame_subset = d['read_only_frame_subset']
    else:
        read_only_frame_subset = False
        

    #For some videos, e.g. some 100x videos, I forgot to get a BG illumination heterogeneity image. For these images, I skip adding the file path to the background image.    
    subtract_bg = True
    if 'subtract_bg' in d_vid:
        if d_vid['subtract_bg'] == 'False' or d_vid['subtract_bg'] == False:
            subtract_bg = False
            print('Skipping background subtraction as directed from the experiment spreadsheet')

    #Load metadata from original image
    v = nd2_handling2.Video_waves(file_path, read_img_directly=False, dir_exp=d_vid['dir_exp'], subtract_bg=subtract_bg, mode=settings_general['mode'])
    v.crop = ap.extract_crop_values_From_df(df_vid) #Crop settings of the raw image
    v.crop2 = ap.extract_crop_values_From_df(df_vid, str='crop_area2')

    if 'frame_to_display' in d:
        frame_to_display = d['frame_to_display']
    else:
        frame_to_display=v.frame_to_display    

    #Set the sub folder which the image stacks are saved:
    if settings_general['hsv_mode'] == 'pixelated':
        sub_folder = 'pixelated' #if super pixel mode is on
    else:
        sub_folder = ''
    print('sub_folder:', sub_folder)

    print('yyy. ',v.n_frames, settings_general['frame_range'])
    #If the set frame range goes beyond the number of frames in the file, set the end of the frame range to the number of frames
    if (settings_general['frame_range'][-1] > v.n_frames) or (settings_general['frame_range'][-1] == 0 and v.n_frames > n_frames_limit and read_only_frame_subset == False):
        settings_general['frame_range'] = range(0,v.n_frames)
        print('settings_general[\'frame_range\']',settings_general['frame_range'])
    v.frame_range = settings_general['frame_range']

    #If the number of frames selected for analysis is larger than the limit (500 frames), split the file up into n_subimgs.
    if len(settings_general['frame_range']) > n_frames_limit:
        split_files_calc_aniso_and_show_output(v,settings_general, d_vid, n_frames_limit, sub_folder, d, p_range = p_range,I_tot_percentiles=I_tot_percentiles, I_tot_limits=I_tot_limits)
    else:
        settings_general['hsv_split_file'] = False
        calc_aniso_and_show_output(v, settings_general, d_vid, d, frame=frame_to_display, sub_folder=sub_folder, p_range = p_range, I_tot_percentiles=I_tot_percentiles, I_tot_limits=I_tot_limits, average_over_frames=average_over_frames)

    #Return an empty dictionary
    dict_var_fun_specific = {}
    return dict_var_fun_specific